Ambulatory or residential? a multi-state analysis of treatments for substance use disorders (Step 1.5)
May 23, 2022
Parametric
Assessment of Fit of Transitions
We need to derive appropriate functional forms and define respective survival functions. One reason to favor patametric models is that they can be easier to generalize. Several candidate distributions were considered to model transitions, including Weibull (assume a monotonically increasing or decreasing hazard), Log-logistic (non-monotonic hazards), Gompertz (assume a monotonically increasing or decreasing hazard), Log-normal (non-monotonic hazards), Exponential (constant hazard), Gamma, Generalized gamma & Generalized F (more flexible).
The following plots fitted survival curves from each model (colored lines), with the Kaplan-Meier estimate, in red, obtained from an example of Jackson available here, added to the contributions of Wathers & Cutler available here.
#options(warn=-1)
n_iter<-10000
tiempo_antes_fits<-Sys.time()
#Weathers, Brandon and Cutler, Richard Dr., "Comparison of Survival Curves Between Cox Proportional
#Hazards, Random Forests, and Conditional Inference Forests in Survival Analysis" (2017). All Graduate
#Plan B and other Reports. 927.
#https://digitalcommons.usu.edu/gradreports/927
#<div style="border: 1px solid #ddd; padding: 5px; overflow-y: scroll; height:350px; overflow-x: scroll; width:100%">
#https://devinincerti.com/2019/01/01/sim-mstate.html
n_trans <- max(trans_matrix, na.rm = TRUE)
fits_wei <- vector(mode = "list", length = n_trans)
fits_weiph <- vector(mode = "list", length = n_trans)
fits_llogis <- vector(mode = "list", length = n_trans)
fits_gomp <- vector(mode = "list", length = n_trans)
fits_logn <- vector(mode = "list", length = n_trans)
fits_exp <- vector(mode = "list", length = n_trans)
fits_gam <- vector(mode = "list", length = n_trans)
fits_ggam <- vector(mode = "list", length = n_trans)
fits_genf <- vector(mode = "list", length = n_trans)
fits_genf_orig <- vector(mode = "list", length = n_trans)
fits_ggam_orig <- vector(mode = "list", length = n_trans)
fits_rp02 <- vector(mode = "list", length = n_trans)
fits_rp03 <- vector(mode = "list", length = n_trans)
fits_rp04 <- vector(mode = "list", length = n_trans)
fits_rp05 <- vector(mode = "list", length = n_trans)
fits_rp06 <- vector(mode = "list", length = n_trans)
fits_rp07 <- vector(mode = "list", length = n_trans)
fits_rp08 <- vector(mode = "list", length = n_trans)
fits_rp09 <- vector(mode = "list", length = n_trans)
fits_rp10 <- vector(mode = "list", length = n_trans)
fits_c_wei <- vector(mode = "list", length = n_trans)
fits_c_weiph <- vector(mode = "list", length = n_trans)
fits_c_llogis <- vector(mode = "list", length = n_trans)
fits_c_gomp <- vector(mode = "list", length = n_trans)
fits_c_logn <- vector(mode = "list", length = n_trans)
fits_c_exp <- vector(mode = "list", length = n_trans)
fits_c_gam <- vector(mode = "list", length = n_trans)
fits_c_ggam <- vector(mode = "list", length = n_trans)
fits_c_genf <- vector(mode = "list", length = n_trans)
fits_c_genf_orig <- vector(mode = "list", length = n_trans)
fits_c_ggam_orig <- vector(mode = "list", length = n_trans)
fits_c_rp02 <- vector(mode = "list", length = n_trans)
fits_c_rp03 <- vector(mode = "list", length = n_trans)
fits_c_rp04 <- vector(mode = "list", length = n_trans)
fits_c_rp05 <- vector(mode = "list", length = n_trans)
fits_c_rp06 <- vector(mode = "list", length = n_trans)
fits_c_rp07 <- vector(mode = "list", length = n_trans)
fits_c_rp08 <- vector(mode = "list", length = n_trans)
fits_c_rp09 <- vector(mode = "list", length = n_trans)
fits_c_rp10 <- vector(mode = "list", length = n_trans)
km.lc <-vector(mode = "list", length = n_trans)
fits_cox<-list()
fits_c_cox<-list()
#"gengamma" Generalized gamma (stable parameterisation)
#"gengamma.orig" Generalized gamma (original parameterisation)
#"genf" Generalized F (stable parameterisation)
#"genf.orig" Generalized F (original parameterisation)
#"weibull" Weibull
#"gamma" Gamma
#"exp" Exponential
#"lnorm" Log-normal
#"gompertz" Gompertz
library(flexsurv)
#Specify the formula
fitform <- Surv(time, status) ~ 1
for (i in 1:n_trans){
fits_wei[[i]] <- flexsurvreg(formula=fitform,
data = subset(ms_d_match_surv, trans == i),
dist = "weibull")
}
for (i in 1:n_trans){
fits_weiph[[i]] <- flexsurvreg(formula=fitform,
data = subset(ms_d_match_surv, trans == i),
dist = "weibullph")
}
for (i in 1:n_trans){
fits_llogis[[i]] <- flexsurvreg(formula=fitform,
data = subset(ms_d_match_surv, trans == i),
dist = "llogis")
}
for (i in 1:n_trans){
fits_gam[[i]] <- flexsurvreg(formula=fitform,
data = subset(ms_d_match_surv, trans == i),
dist = "gamma")
}
#LOS ERRORES OCURREN CUANDO NO HAY COVARIABLES
#In (function (q, shape, rate = 1, scale = 1/rate, lower.tail = TRUE, ... : NaNs produced
#gamma no funcionó: NaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs produced
for (i in 1:n_trans){
fits_ggam[[i]] <- flexsurvreg(formula=fitform,
data = subset(ms_d_match_surv, trans == i),
dist = "gengamma")
}
for (i in 1:n_trans){
fits_gomp[[i]] <- flexsurvreg(formula=fitform,
data = subset(ms_d_match_surv, trans == i),
dist = "gompertz")
}## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.32590883280109: valor inicial en 'vmmin' no es finitofor (i in 1:n_trans){
fits_logn[[i]] <- flexsurvreg(formula=fitform,
data = subset(ms_d_match_surv, trans == i),
dist = "lnorm")
}
for (i in 1:n_trans){
fits_exp[[i]] <- flexsurvreg(formula=fitform,
data = subset(ms_d_match_surv, trans == i),
dist = "exp")
}
no_attempts <- 20
for (i in 1:n_trans){
r <- NULL
attempt <- 0
while( is.null(r) && attempt <= no_attempts ) {
attempt <- attempt + 1
try(
r <- flexsurvreg(formula=fitform,
data = subset(ms_d_match_surv, trans == i),
dist = "gengamma.orig")
)
}
fits_ggam_orig[[i]] <- r
}
for (i in 1:n_trans){
r <- NULL
attempt <- 0
while( is.null(r) && attempt <= no_attempts ) {
attempt <- attempt + 1
try(
r <- flexsurvreg(formula=fitform,
data = subset(ms_d_match_surv, trans == i),
dist = "genf.orig")
)
}
fits_genf_orig[[i]] <- r
} ## Warning in flexsurvreg(formula = fitform, data = subset(ms_d_match_surv, :
## Optimisation has probably not converged to the maximum likelihood - Hessian is
## not positive definite.for (i in 1:n_trans){
r <- NULL
attempt <- 0
while( is.null(r) && attempt <= no_attempts ) {
attempt <- attempt + 1
try(
r <- flexsurvreg(formula=fitform,
data = subset(ms_d_match_surv, trans == i),
dist = "genf")
)
}
fits_genf[[i]] <- r
}
for (i in 1:n_trans){
r <- NULL
attempt <- 0
while( is.null(r) && attempt <= no_attempts ) {
attempt <- attempt + 1
try(
r <- flexsurvspline(formula=fitform,k=1,
data = subset(ms_d_match_surv, trans == i))
)
}
fits_rp02[[i]] <- r
}
for (i in 1:n_trans){
r <- NULL
attempt <- 0
while( is.null(r) && attempt <= no_attempts ) {
attempt <- attempt + 1
try(
r <- flexsurvspline(formula=fitform,k=2,
data = subset(ms_d_match_surv, trans == i))
)
}
fits_rp03[[i]] <- r
}
for (i in 1:n_trans){
r <- NULL
attempt <- 0
while( is.null(r) && attempt <= no_attempts ) {
attempt <- attempt + 1
try(
r <- flexsurvspline(formula=fitform,k=3,
data = subset(ms_d_match_surv, trans == i))
)
}
fits_rp04[[i]] <- r
}
for (i in 1:n_trans){
r <- NULL
attempt <- 0
while( is.null(r) && attempt <= no_attempts ) {
attempt <- attempt + 1
try(
r <- flexsurvspline(formula=fitform,k=4,
data = subset(ms_d_match_surv, trans == i))
)
}
fits_rp05[[i]] <- r
}
for (i in 1:n_trans){
r <- NULL
attempt <- 0
while( is.null(r) && attempt <= no_attempts ) {
attempt <- attempt + 1
try(
r <- flexsurvspline(formula=fitform,k=5,
data = subset(ms_d_match_surv, trans == i))
)
}
fits_rp06[[i]] <- r
}
for (i in 1:n_trans){
r <- NULL
attempt <- 0
while( is.null(r) && attempt <= no_attempts ) {
attempt <- attempt + 1
try(
r <- flexsurvspline(formula=fitform,k=6,
data = subset(ms_d_match_surv, trans == i))
)
}
fits_rp07[[i]] <- r
}
for (i in 1:n_trans){
r <- NULL
attempt <- 0
while( is.null(r) && attempt <= no_attempts ) {
attempt <- attempt + 1
try(
r <- flexsurvspline(formula=fitform,k=7,
data = subset(ms_d_match_surv, trans == i))
)
}
fits_rp08[[i]] <- r
}
for (i in 1:n_trans){
r <- NULL
attempt <- 0
while( is.null(r) && attempt <= no_attempts ) {
attempt <- attempt + 1
try(
r <- flexsurvspline(formula=fitform,k=8,
data = subset(ms_d_match_surv, trans == i))
)
}
fits_rp09[[i]] <- r
}
for (i in 1:n_trans){
r <- NULL
attempt <- 0
while( is.null(r) && attempt <= no_attempts ) {
attempt <- attempt + 1
try(
r <- flexsurvspline(formula=fitform,k=9,
data = subset(ms_d_match_surv, trans == i))
)
}
fits_rp10[[i]] <- r
}
for (i in 1:n_trans){
km.lc[[i]] <- survfit(formula= fitform, data = subset(ms_d_match_surv, trans == i))
}
transition_label<- plots$title
attr(transition_label,"names")<- plots$trans
get_distinct_hues <- function(ncolor,s=0.5,v=0.95,seed=2125) {
golden_ratio_conjugate <- 0.618033988749895
set.seed(seed)
h <- runif(1)
H <- vector("numeric",ncolor)
for(i in seq_len(ncolor)) {
h <- (h + golden_ratio_conjugate) %% 1
H[i] <- h
}
hsv(H,s=s,v=v)
}distinct_hues<-get_distinct_hues(20)
layout(matrix(1:n_trans, nc = 2, byrow = F))
for (i in 1:n_trans){
plot(km.lc[[i]], col="red", conf.int=F);
lines(fits_wei[[i]], col=distinct_hues[1], ci=F);
lines(fits_weiph[[i]], col=distinct_hues[2], ci=F);
lines(fits_gomp[[i]], col=distinct_hues[3], ci=F);
lines(fits_llogis[[i]], col=distinct_hues[4], ci=F);#A0A36D
lines(fits_gam[[i]], col=distinct_hues[5], ci=F);
lines(fits_ggam[[i]], col=distinct_hues[6], ci=F);
lines(fits_ggam_orig[[i]], col=distinct_hues[7], ci=F);
lines(fits_logn[[i]], col=distinct_hues[8], ci=F);
lines(fits_exp[[i]],col=distinct_hues[9], ci=F);
lines(fits_genf[[i]],col=distinct_hues[10], ci=F, lty = 2);
lines(fits_genf_orig[[i]],col=distinct_hues[11], ci=F, lty = 2);
lines(fits_rp02[[i]],col=distinct_hues[12], ci=F, lty = 2);
lines(fits_rp03[[i]],col=distinct_hues[13], ci=F, lty = 2);
lines(fits_rp04[[i]],col=distinct_hues[14], ci=F, lty = 2);
lines(fits_rp05[[i]],col=distinct_hues[15], ci=F, lty = 2);
lines(fits_rp06[[i]],col=distinct_hues[16], ci=F, lty = 2);
lines(fits_rp07[[i]],col=distinct_hues[17], ci=F, lty = 3);
lines(fits_rp08[[i]],col=distinct_hues[18], ci=F, lty = 3);
lines(fits_rp09[[i]],col=distinct_hues[19], ci=F, lty = 3);
lines(fits_rp10[[i]],col=distinct_hues[20], ci=F, lty = 3)
legend("bottomleft", legend = c("Kaplan-Meier","Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
"Generalized gamma", "Generalized gamma (orig)", "Lognormal", "Exponential", "Generalized F","Generalized F (orig)","RP2","RP3","RP4","RP5", "RP6","RP7", "RP8", "RP9", "RP10"), col =
c("red",distinct_hues),
title = "Distributions", cex = .95, bty = "n", lty=c(rep(1,10),rep(2,6),rep(3,4)),lwd=3)# lty = 1:2,
title(main=transition_label[[i]])
}
Figure 19. Vissual Assessment of Survival Curves
endTime <- Sys.time()
paste0("Time in process: ");endTime - tiempo_antes_fits## [1] "Time in process: "## Time difference of 3.497949 mins#23 min aprox.
#For more complicated models, users should specify what covariate values they want summaries for, rather than relying on the default
#</div>
options(warn=0)
if(no_mostrar==1){
jpeg(paste0(dta_path,"_mult_state_ags/exp_surv_int_only.jpg"), height=30, width= 30, res= 600, units = "in")
layout(matrix(1:n_trans, nc = 2, byrow = F))
for (i in 1:n_trans){
plot(km.lc[[i]], col="red", conf.int=F);
lines(fits_wei[[i]], col=distinct_hues[1], ci=F);
lines(fits_weiph[[i]], col=distinct_hues[2], ci=F);
lines(fits_gomp[[i]], col=distinct_hues[3], ci=F);
lines(fits_llogis[[i]], col=distinct_hues[4], ci=F);#A0A36D
lines(fits_gam[[i]], col=distinct_hues[5], ci=F);
lines(fits_ggam[[i]], col=distinct_hues[6], ci=F);
lines(fits_ggam_orig[[i]], col=distinct_hues[7], ci=F);
lines(fits_logn[[i]], col=distinct_hues[8], ci=F);
lines(fits_exp[[i]],col=distinct_hues[9], ci=F);
lines(fits_genf[[i]],col=distinct_hues[10], ci=F, lty = 2);
lines(fits_genf_orig[[i]],col=distinct_hues[11], ci=F, lty = 2);
lines(fits_rp02[[i]],col=distinct_hues[12], ci=F, lty = 2);
lines(fits_rp03[[i]],col=distinct_hues[13], ci=F, lty = 2);
lines(fits_rp04[[i]],col=distinct_hues[14], ci=F, lty = 2);
lines(fits_rp05[[i]],col=distinct_hues[15], ci=F, lty = 2);
lines(fits_rp06[[i]],col=distinct_hues[16], ci=F, lty = 2);
lines(fits_rp07[[i]],col=distinct_hues[17], ci=F, lty = 3);
lines(fits_rp08[[i]],col=distinct_hues[18], ci=F, lty = 3);
lines(fits_rp09[[i]],col=distinct_hues[19], ci=F, lty = 3);
lines(fits_rp10[[i]],col=distinct_hues[20], ci=F, lty = 3)
legend("bottomleft", legend = c("Kaplan-Meier","Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
"Generalized gamma", "Generalized gamma (orig)", "Lognormal", "Exponential", "Generalized F","Generalized F (orig)","RP2","RP3","RP4","RP5", "RP6","RP7", "RP8", "RP9", "RP10"), col =
c("red",distinct_hues),
title = "Distributions", cex = .95, bty = "n", lty=c(rep(1,10),rep(2,6),rep(3,4)),lwd=3)# lty = 1:2,
title(main=transition_label[[i]])
}
dev.off()
}#
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
fit_flexsurvreg0<-data.frame()
fitted_flexsurvreg0<-data.frame()
dists_no_covs_10s<-cbind.data.frame(covs=c(rep("fits_",(20)*n_trans)),
formal=rep(c("Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
"Generalized gamma", "Generalized gamma (original)", "Lognormal", "Exponential", "Generalized F", "Generalized F(original)", paste0("RP0",2:9),"RP10"),1*n_trans),
dist=c("weibull", "weibullph", "gompertz", "llogis", "gamma", "gengamma","gengamma.orig", "lnorm", "exp", "genf","genf.orig",rep("no dist",9)),
model=rep(c("wei", "weiph", "gomp", "llogis", "gam","ggam", "ggam_orig", "logn", "exp", "genf","genf_orig", paste0("rp0",2:9),"rp10"),1*n_trans),
trans=rep(1:n_trans, each=20))
for (i in 1:nrow(dists_no_covs_10s)){ #
cat(paste0("#### Flexible Survival Model (w/ covs): ",
dists_no_covs_10s[i,"formal"], "; transition: ",dists_no_covs_10s[i,"trans"], "\n \n"))
model<-paste0("fits_",dists_no_covs_10s[i,"model"])
if(!is.null(get(model)[[dists_no_covs_10s[i,"trans"]]])){
fitted_flexsurvreg0<- rbind(fitted_flexsurvreg0,cbind.data.frame(dist=rep(dists_no_covs_10s[i,"formal"],),
trans=rep(dists_no_covs_10s[i,"trans"],),
data.table::data.table(summary(get(model)[[dists_no_covs_10s[i,"trans"]]],
#newdata= newdat2a,
#newtime=unique(fitted_km$time),
type = "survival",
tidy=T))))
# Generate fit indices
fit_flexsurvreg0<-rbind(fit_flexsurvreg0,
cbind(dist= dists_no_covs_10s[i,"formal"],
transition=dists_no_covs_10s[i,"trans"],
fitstats.flexsurvreg(get(model)[[dists_no_covs_10s[i,"trans"]]])))
#the BIC may not be appropriate if none of the candidate models are considered to be close to the ‘true’ model.
} else {
cat(paste0("The model that assumed a ",dists_no_covs_10s[i,"formal"]," distribution for the transition number ",dists_no_covs_10s[i,"trans"]," did not converge \n\n"))
}
}## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 1
##
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 1
##
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 1
##
## The model that assumed a Gompertz distribution for the transition number 1 did not converge
##
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 1
##
## #### Flexible Survival Model (w/ covs): Gamma; transition: 1
##
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 1
##
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 1
##
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 1
##
## #### Flexible Survival Model (w/ covs): Exponential; transition: 1
##
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 1
##
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 1
##
## #### Flexible Survival Model (w/ covs): RP02; transition: 1
##
## #### Flexible Survival Model (w/ covs): RP03; transition: 1
##
## #### Flexible Survival Model (w/ covs): RP04; transition: 1
##
## #### Flexible Survival Model (w/ covs): RP05; transition: 1
##
## #### Flexible Survival Model (w/ covs): RP06; transition: 1
##
## #### Flexible Survival Model (w/ covs): RP07; transition: 1
##
## #### Flexible Survival Model (w/ covs): RP08; transition: 1
##
## #### Flexible Survival Model (w/ covs): RP09; transition: 1
##
## #### Flexible Survival Model (w/ covs): RP10; transition: 1
##
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 2
##
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 2
##
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 2
##
## The model that assumed a Gompertz distribution for the transition number 2 did not converge
##
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 2
##
## #### Flexible Survival Model (w/ covs): Gamma; transition: 2
##
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 2
##
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 2
##
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 2
##
## #### Flexible Survival Model (w/ covs): Exponential; transition: 2
##
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 2
##
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 2
##
## #### Flexible Survival Model (w/ covs): RP02; transition: 2
##
## #### Flexible Survival Model (w/ covs): RP03; transition: 2
##
## #### Flexible Survival Model (w/ covs): RP04; transition: 2
##
## #### Flexible Survival Model (w/ covs): RP05; transition: 2
##
## #### Flexible Survival Model (w/ covs): RP06; transition: 2
##
## #### Flexible Survival Model (w/ covs): RP07; transition: 2
##
## #### Flexible Survival Model (w/ covs): RP08; transition: 2
##
## #### Flexible Survival Model (w/ covs): RP09; transition: 2
##
## #### Flexible Survival Model (w/ covs): RP10; transition: 2
##
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 3
##
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 3
##
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 3
##
## The model that assumed a Gompertz distribution for the transition number 3 did not converge
##
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 3
##
## #### Flexible Survival Model (w/ covs): Gamma; transition: 3
##
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 3
##
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 3
##
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 3
##
## #### Flexible Survival Model (w/ covs): Exponential; transition: 3
##
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 3
##
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 3
##
## #### Flexible Survival Model (w/ covs): RP02; transition: 3
##
## #### Flexible Survival Model (w/ covs): RP03; transition: 3
##
## #### Flexible Survival Model (w/ covs): RP04; transition: 3
##
## #### Flexible Survival Model (w/ covs): RP05; transition: 3
##
## #### Flexible Survival Model (w/ covs): RP06; transition: 3
##
## #### Flexible Survival Model (w/ covs): RP07; transition: 3
##
## #### Flexible Survival Model (w/ covs): RP08; transition: 3
##
## #### Flexible Survival Model (w/ covs): RP09; transition: 3
##
## #### Flexible Survival Model (w/ covs): RP10; transition: 3
##
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 4
##
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 4
##
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 4
##
## The model that assumed a Gompertz distribution for the transition number 4 did not converge
##
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 4
##
## #### Flexible Survival Model (w/ covs): Gamma; transition: 4
##
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 4
##
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 4
##
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 4
##
## #### Flexible Survival Model (w/ covs): Exponential; transition: 4
##
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 4
##
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 4
##
## #### Flexible Survival Model (w/ covs): RP02; transition: 4
##
## #### Flexible Survival Model (w/ covs): RP03; transition: 4
##
## #### Flexible Survival Model (w/ covs): RP04; transition: 4
##
## #### Flexible Survival Model (w/ covs): RP05; transition: 4
##
## #### Flexible Survival Model (w/ covs): RP06; transition: 4
##
## #### Flexible Survival Model (w/ covs): RP07; transition: 4
##
## #### Flexible Survival Model (w/ covs): RP08; transition: 4
##
## #### Flexible Survival Model (w/ covs): RP09; transition: 4
##
## #### Flexible Survival Model (w/ covs): RP10; transition: 4
## #_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
if(no_mostrar==1){
fit_flexsurvreg0 %>%
dplyr::group_by(transition) %>%
summarise(mean(ucl,na.rm=T))
}
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
km.lc0<-list()
fitted_km0<-data.frame()
for (i in 1:n_trans){
km.lc0[[i]] <- survfit(formula= fitform, data = subset(ms_d_match_surv, trans == i))
fitted_km0<-rbind(fitted_km0,cbind.data.frame(trans=i,fortify(km.lc0[[i]])))
}
#Calculate error
#c("Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma","Generalized gamma", "Lognormal", "Exponential", "Generalized F")
fitted_flexsurvreg_binned_mix0<-
data.frame(fitted_flexsurvreg0)[,c("dist","trans","time","est","lcl","ucl")] %>%
dplyr::left_join(fitted_km0, by=c("trans","time")) %>%
#there are many observations that was not available due to empirical missing data
#dplyr::filter(!is.na(surv))
dplyr::group_by(dist,trans) %>%
#dplyr::mutate(est= ifelse(is.na(est), mean(est, na.rm=TRUE), est)) %>%
#dplyr::mutate(surv= ifelse(is.na(surv), mean(surv, na.rm=TRUE), surv)) %>%
tidyr::fill(surv,.direction="updown") %>%
tidyr::fill(est,.direction="updown") %>%
ungroup()
db_for_apply_rmse0<-
cbind.data.frame(
dist=rep(c("Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
"Generalized gamma", "Generalized gamma (original)", "Lognormal", "Exponential", "Generalized F", "Generalized F (orignial)", paste0("RP0",2:9),"RP10"),n_trans),
trans=rep(c(1:n_trans),each=20*2))
rmse_comp_fits0<- data.frame()
for(i in 1:nrow(db_for_apply_rmse0)){
rmse<- Metrics::rmse(subset(fitted_flexsurvreg_binned_mix0, dist==db_for_apply_rmse0[i,"dist"] &
trans==db_for_apply_rmse0[i,"trans"])$est,
subset(fitted_flexsurvreg_binned_mix0, dist==db_for_apply_rmse0[i,"dist"] &
trans==db_for_apply_rmse0[i,"trans"])$surv)
rmse_comp_fits0<- rbind(rmse_comp_fits0,cbind(dist=db_for_apply_rmse0[i,"dist"],
residential=db_for_apply_rmse0[i,"residential"],
trans=db_for_apply_rmse0[i,"trans"],
rmse=rmse))
}
rmse_comp_fits0_filter<-
rmse_comp_fits0 %>%
dplyr::filter(rmse!="NaN") %>%
dplyr::arrange(trans,dist, rmse)%>%
dplyr::mutate(rmse=round(as.numeric(rmse),4))
setting_type_label<- c(`0`="Outpatient",`1`="Residential")
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_We checked the AICs, cumulative relative error (differences between the Kaplan-Meier curve, and the observed behavior of each distribution).
c26 <- c(
"dodgerblue2", # red
"green4",
"#6A3D9A", # purple
"#FF7F00", # orange
"gray16", "gold1",
"skyblue2", "#FB9A99", # lt pink
"palegreen2",
"#CAB2D6", # lt purple
"#FDBF6F", # lt orange
"gray70", "khaki2",
"maroon", "orchid1", "deeppink1", "blue1", "steelblue4",
"darkturquoise", "green1", "yellow4", "yellow3",
"darkorange4", "brown", "gray40")
fitted_flexsurvreg_binned_mix0_plot<-
fitted_flexsurvreg_binned_mix0 %>%
dplyr::mutate(abs((est-surv)/surv)) %>% group_by(dist,trans) %>% dplyr::mutate(cumsum_error=cumsum(`abs((est - surv)/surv)`)) %>%
dplyr::ungroup() %>%
dplyr::mutate(text=paste0("Distribution: ",dist,"<br>survival: ",sprintf("%1.2f",est),"<br>Time (in years): ",round(time/365.25,2), "<br>CIs: ",ifelse(!is.na(lcl),"CIs","no CIs"),"<br>Cumsum(abs((est-surv)/surv))= ",round(cumsum_error,2))) %>%
dplyr::mutate(text2=paste0("Distribution: KM","<br>survival: ",sprintf("%1.2f",surv),"<br>Time (in years): ",round(time/365.25,2))) %>%
ggplot()+
geom_line(aes(time, est, color=dist),size=1)+
geom_line(aes(time, surv), color="red",size=1)+
scale_color_manual(name="Distributions", values = c26[1:20])+
#c("black",brewer.pal(n = 9, name = 'Paired')))+
#c("#112A60","#085754","#D3A347","#4F3C91","red","#112A60","#085754","#8F630D","#251363")) +
facet_wrap(.~trans,labeller = labeller(trans = transition_label), ncol=3, scales = "free_y")+
theme_bw()+
theme(legend.position="bottom",
strip.background = element_rect(fill = "white", colour = "white"))+
scale_x_continuous(breaks = seq(0, 365.25*12, 2*365.25), labels=seq(0,12,2))+
#ylim(0,1.25)+
#theme(axis.text.x = element_blank(),
# panel.grid.major = element_blank(),
# panel.grid.minor = element_blank()) +
labs(y="",x="")+
theme(legend.position='none')
m <- list(
l = 80*1.05,
r = 80,
b = 80,
t = 80*1.05,
pad = 4
)
p1<-
fitted_flexsurvreg_binned_mix0_plot$data %>%
dplyr::mutate(years=round(time/365.25,0)) %>%
dplyr::mutate(trans2=trans) %>%
dplyr::left_join(cbind.data.frame(transition_label,trans2=1:n_trans), by="trans2") %>%
group_by(trans) %>%
group_map(~ plot_ly(data=., x = ~time, y = ~est, color = ~dist, type = "scatter", mode="lines", text=~text) %>%
add_trace(y = ~surv, type = 'scatter', mode = 'lines', line = list(color = 'red', width = 1), text=~text2)%>%
layout(annotations = list(list(x = 0.5 , y = 1.03, text = ~unique(transition_label), showarrow = F, xref='paper', yref='paper'))) %>%
layout(
xaxis = list(
ticktext = list(1, 3, 5, 7, 9, 11),
textangle = 0,
tickvals = list(365.25*1, 365.25*3, 365.25*5, 365.25*7, 365.25*9, 365.25*11),
tickmode = "array"
)), keep=TRUE) %>%
subplot(nrows = 3, shareX = T, shareY=T, titleX = F, titleY = F, margin = .05)%>% layout(showlegend = F) %>% #, margin = 0.05) %>%
layout(annotations = list(
list(x = -0.03, y = 0.5, text = "Survival",
font = list(color = "black",size = 15),
textangle = 270,
showarrow = F, xref='paper', yref='paper', size=48,margin =m)))
p1